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p-"] Abstract 



Under the effect of strong genetic drift, it is highly probable to observe gene 
fixation or gene loss in a population, shown by infinite peaks on a coherently 
constructed potential energy landscape. It is then important to ask what such 
singular peaks imply, with or without the effects of other biological factors. 
We studied the stochastic escape time from the infinite potential peaks in 
the Wright-Fisher model, where the typical two-scale diffusion dynamics was 
observed via computer simulations. We numerically found the average escape 
time for all the bi-stable cases and analytically approximated the results 
under weak mutations and selections by calculating the mean first passage 
time (MFPT) in singular potential peak. Our results showed that Kramers' 
classical escape formula can be extended to the models with non-Gaussian 
probability distributions, overcoming constraints in previous methods. The 
constructed landscape provides a global and coherent description for system's 
;^j evolutionary dynamics, allowing new biological results to be generated. 
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1. Introduction 



Interactions between different biological driving forces can generate very 
complex phenomena in evolution. These forces (genetic drift, mutation, se- 
lection, etc.) vary in their forms and intensities of effects, characterized by 



their different operating timescales (Gillespie, 1998). In the study of biologi 



cal evolution, one of the most important and interesting issues is to describe 



and separate the multiple timescale dynamics (Gillespie, 1984) and to an 



alyze the occurring rates of rare events (Kimura, 1957). Related problems 



have been referred in different contexts in literature, in close relation with the 



concept of adaptive landscape (Wright, 1932 Arnold et al., 2001; Ao, 2005). 



In population genetics, adaptation was found not to be limited by the rate 



toward local adaptive peaks but by the peak-to-peak transition rate (Wright 



1932). Studies were carried on how the success or failure of a mutant gene 



depends on chance for all levels of selective dominance (Kimura, 1962). A 



knowledge of the frequencies with which populations move toward a local 
peak and that move between different peaks helps understand the mecha- 



nisms of adaptation and divergence (Barton and Rouhani, 1987). It was also 



reviewed how genetic barriers for gene flow are established and related it to 
biological speciation (Gavrilets, 2003). Results on multiple adaptive peaks 



and the associated multiple evolutionary timescales were found important for 



studying evolutionary robustness (Ao, 2009). The ideas and methodologies 



are also widely discussed outside biology (Qian [2005 ; Kriiger, 2010[ ) 

Typically, the existence of a genetic barrier (or adaptive valley) suggests 
the separation of different evolutionary timescales. In chemistry, the known 
Arrhenius formula estimates the separation factor to be an exponential term 
of the energy barrier height (or valley depth) (Hanggi et al. , 1990). It was 



latter systematically studied in thermally activated systems (Kramers, 1940) 



In population genetics, however, random drift may cause problems for the 
the classical results. The key point here is that the strength of drift-induced 
noise varies from state to state and vanishes at the fixation states. This 
is different from the usual assumption of constant weak noise in physical 



processes (Gardiner, 1985), and the system probability distribution is often 



far from Gaussian near local equilibria. 

To investigate the effects of genetic drift and its interactions with other 
biological factors, we study the Wright-Fisher diffusion process, a classical 



model for the study of genetic drift and multi-scale dynamics (Kimura, 1964 



Blythe and McKane, 2007). It assumes fixed population size and considers 



the change of portions of different gene copies. In this model, we construct 
a potential function which is exactly consistent to the potential energy in 
physical sciences. It can be visualized as a potential landscape, in compar- 
ison to the classical fitness landscape (Wright, 1932). Under strong genetic 
drift, a population is very likely to be found near a potential peak at the 
fixation state. The tendency is shown by singular (infinite) peaks on the 
potential landscape at such states. The questions here are: What is the bio- 
logical meaning of these infinite peaks? Do they imply the ultimate fixation 
of a gene type? If not, what is the life time (average escape time) of such 
states? In such cases, the classical escape formula gives biologically unex- 
pected estimations for the results. The escape times are often estimated by 
calculating the numerical solutions of the mean first passage times (MFPT), 



but the relation between the two concepts was not made clear (Kimura and 



Ohta. 1969 



Karlin and Taylor, 1981 Lande, 1985). In this article, we try 



to formally find the relation between the escape time and the MFPT in the 
present model. We look for the analytical estimation of the escape time from 
the singular potential peaks, and see whether such peaks imply gene fixation. 
In more physical contexts, the stochastic escape time in a diffusion pro- 
cess were studied by different methods. One such example is to calculate the 



stationary probability flux from an attractive basin (Kramers, 1940), known 



as the flux-over-population method. In general diffusion process, this method 
is shown equivalent to the MFPT calculations (Hanggi et al. , 1990). It often 



assumes Gaussian-like probability distribution near an adaptive peak, which 
is not the fact in the present model. Another approach is to calculate the 
eigenvalue of the diffusion equations (Barton and Rouhani, 1987). However, 



the method is not always applicable for bi-stable population models, espe- 
cially when selection is weak. A third approach is more from the side of 
population genetics, known as the "rate of genetic substitution" (Kimura 



1962 Gillespie, 1998). It bypasses the real dynamics in infinite potential 



peaks, but the results are restricted to the models where fixation is possible 
and mutation is very weak. In the present article, we study the simulation 
results of the discrete Wright-Fisher model and approximate the average es- 
cape time by calculating the MFPT. Based on previous results (Xu et al 



2012), we study the relation of escape time and MFPT under the present 
framework. Our results show that Kramers' classical escape formula can be 
extended to the models with singular potential landscape, overcoming previ- 
ous constraints and allowing more complex dynamics. Our results provide a 
complete answer for the bi-stable problems in the present model. 



The present article is organized as follows: In Section 2, we introduce 
the 1-d diffusion process and define a potential landscape. In Section 3, 
we discuss the uphill and downhill landscape dynamics under strong genetic 
drift. In Section 4, we first simulate the discrete Wright-Fisher model and 
analyze the simulated escape rate. We then analytically approximate the 
average escape time by calculating the MFPT in the infinite potential. After 
that, we come to two models with selections. In Section 5, we discuss the 
relation between the MFPT and the escape time. We then compare our 
results of the escape time and our potential landscape to others' work in 
literature. To conclude, we discuss the fixation condition in infinite potential 
and other biological insights based on the present work. 



2. Wright-Fisher model and potential landscape 

2. 1 . Diffusion process 

The 1-d Wright-Fisher model considers the evolution of a diploid popu- 
lation at one locus. The number of individuals in a population is fixed at 
N and the generations are non-overlapping. Denote the interested pair of 
alleles as A\ and A 2 ; the total number of the gene copies of A\ and A 2 in 
the population gene pool is 2N. In the present article, we mainly study the 
continuous diffusion approximation of the Wright-Fisher model (assume N is 
big enough for the continuous approximation). Let the frequency of A\ gene 
be x, so the frequency of A 2 is 1 — x. Let p(x, t) be the probability distribu- 
tion of A\ at time t. The diffusion equation for the continuous Wright-Fisher 



model is given by (Kimura, 1964 Ewens, 2004 Blythe and McKane, 2007): 



d t p{x, t) = -d 2 x 



V(x)p(x,t) 



d x 



M(x)p(x,t) 



(1) 



M (x) is the average change of the A\ frequency per generation, correspond- 
ing to the deterministic factors of the system. V(x) is the variance of the 
stochastic factors. For example, under mutation and selection: 



M(x) = -px + z/(l - x) + 



x(l — x) du 



(2) 



2w dx 

where p is the mutation rate from A\ to A 2 , v is that from A 2 X>o A\] u gives 
the average fitness of the population in that generation, which depends on 
x. Under random genetic drift: 



V{x) =x{l-x)/2N . 



(3) 



The population size 2N (the number of gene copies on the considered locus 
in the diploid population gene pool) controls the intensity of genetic drift. 

In 1-d model, by assuming zero probability current at system equilibrium 
(t = +00), the equilibrium probability distribution of Eq. ((!]) can be easily 
obtained as (Gardiner, 1985): 



p(x, +00) 



1 



V{x) 



cxp 



2M(y) 

V(y) 



dy 



cxp 



2M(y) - V'(y) 

V(y) 



dy 



where the normalization constant Z is given by 



exp 



2M(y) - V'(y) 

V(y) 



dy 



dx . 



(4) 



(5) 



2.2. Potential landscape 

The form of Eq. Q immediately suggests the existence of a potential 
function <&{x) from the Boltzmann-Gibbs distribution in physics (Ao, 2005): 



p(x, +00) oc exp($(x)) . 



(6) 



This definition directly connects the equilibrium probability distribution to a 
potential energy function, which can be visualized as a potential landscape: 



$(x) 



* 2M(y) - V'{y) 

V(y) 



dy 



x m 

D(y) 



dy . 



(7) 



The maxima and minima states on the equilibrium distribution are exactly 
the peaks and valleys on the potential landscape. Here we have defined a 
directed force f(x) and an undirected diffusion term D(x), which are closely 
related to the system's long-term dynamics: 



f(x) = M(x) - V'(x)/2 
D(x) = V(x)/2 . 



(8) 
(9) 



We can specify Eq. ([7]) in the Wright-Fisher model under the effects of 
genetic drift, mutation, and selection, by considering Eqs. ^ and ^: 



$(x) 



lnx(l -x) +4N 



vhax + /iln(l 



x) 



+ 2iVlncJ. 



(10) 



With the analytical form of $(x), we may classify the Wright-Fisher diffusion 
models under different parameters according to their long-term behaviors. 
Several examples are given in Figure 1. 



3. Bi-stable dynamics under strong genetic drift 

3.1. Mutation and genetic drift 

We apply above methods in the bi-stable cases in the Wright-Fisher 
model. We start with the simplest mutation-drift case, 

$(x) = (4Nv - 1) lnx + (4JV/x - 1) ln(l - x) . (11) 

To maintain a bi-stable system, we set 4iW, 4iV/i < 1. There is a unique 
valley state (saddle point) in (0, 1), here we denote as x = a: 

a= (l-4AV)/(2 -AN/i-ANi/) , (12) 

satisfying $'( a ) = an d <£"( a ) > 0. Such a valley defines two attractive 
basins (0, a) and (a, 1). A population starting in (0, a) (or (a, 1)) is expected 
(averaging the effect of noise) to reach x = (or x — 1) peak as a result 
of directed evolution driven by the directed force /. Compared to the usual 
adaptive force induced by selection, the directed force integrates the effects 
of other factors (here genetic drift and mutation) besides selection (if there 
is). Under mutation and drift, there is 

f(x) = -fix + i/(l - x) - (1 - 2x)/4N . (13) 

Obviously there is f(a) = 0. We will have more detailed discussion on this 
point in the next section. The potential landscapes under different parame- 
ters are shown in Figure 1. 

3.2. Uphill dynamics 

The U-shaped bi-stable landscapes in Figure 1 immediately suggest the 
existence of two-timescale dynamics. The movements of a population, if vi- 
sualized on the potential landscape, can be classified into two fundamentally 
different types: uphill and downhill processes, often demonstrated to oper- 



ate on two different timescales (Zhou and Qian, 2011). In Figure 2, as an 



example, a population starting in the (0, a) basin reaches x = (uphill) in 
71, which is dominated by the characteristic time of the uphill movements; it 
then escapes to x = 1 in 72, which is dominated by the downhill movements 
— > a. The usual assumption is 71 ^C 7^- 

The uphill valley-to-peak evolution in 7i is mainly driven by the directed 
forces, here our f(x). We refer to the Langevin equation that describes the 
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same evolutionary process with Eq. (|T|), but from the point of view of a single 
population's stochastic evolution ( |Ao et aL] 2007): 



x 



f(x) + ^Dix~)C(x,t) 



(14) 



Here x = dx/dt denotes the change rate of x. ({x,t) is the Gaussian white 
noise. Note that Eq. (14) is related to Eq . (|T|) via a different stochastic 
integral from those of Ito and Stratonovich (Ao et al. , 2007). By averaging 



the effects of noise over its probability distribution (instead of taking the 
zero-noise limit A" —¥ +oo), we obtain the average uphill rate: 



x = f(x) . 



(15) 



It is easy to verify that $ is non-decreasing along the noise-free evolutionary 
trajectory of a population: 



$ = $'(x) ■ x = f{x)/D{x) > . 



(16) 



This manifests Wright's essential idea of non-decreasing scalar function that 
can be visualized as a potential landscape. For linear f(x), we can always 
take the approximation form / ~ — \f\x (here we replace x with x — a; note 
that x gives the distance between x and a). The solution of Eq. (15) takes 
the approximate form 



x(t) = x(0) • exp(-|/|£) = x(0) • exp(-t/Ti) 



(17) 



where x(0) gives the initial state of the population, and 7i is usually called 
the relaxation time. Under strong genetic drift (4A/V, 4A 7 // ^C 1), the uphill 
rate is 

f{x) w -(1 - 2x)/AN . (18) 



x 



For x < 1/2, a population is expected to move toward x = 0, and vice versa. 
This is consistent to the biological expectation, that a population is expected 
to be fixed at either monomorphic state (all the individuals have the same 
gene type A\A\, or all are A 2 A 2 ) The first time scale for uphill dynamics to 
reach the local potential peak is 



T 1 ~\f\- 1 = 2N-0(l) 



(19) 



the typical operating timescale of genetic drift (Gillespie, 1998). 



3.3. Downhill dynamics 

The downhill dynamics has not an explicit evolutionary trajectory, but 
is considered as accumulations of rare downhill movements. The process 
can be characterized by the waiting time r for such rare effects to grow big 
enough, so that a population escapes from the original attractive basin, goes 
through the potential valley, and stays stable in another attractive basin. In 
a diffusion model with a finite potential barrier (potential valley), a classical 
formula estimates the escape time as ( |Kramers 1940) 

(20) 



T x exp (A$) . 



Here A<J> is the potential barrier height (valley depth), exp (A<3>) is the 
Arrhenius term. 

For the downhill dynamics in the present case, the classical formula would 
give an estimation of infinite escape time under strong genetic drift and 
weak mutation: From Eq. (11), there is $(0) = +oo, which leads to A$ = 
$(0) - $(a) = +oo and thus by Eq. (J20l) 



r 



-oo 



(21) 



Under pure genetic drift, this infinite escape time is expected. A pop- 
ulation will finally be fixed at either A\ or A 2 gene. Without mutation or 
other input of different genes, the fixation state will not be changed. This 
is also shown by the stationary distribution as a combination of the Dirac 



delta functions under pure genetic drift (McKane and Waxman, 2007): 



p(x, oo) = (1 - C)5(x) + C8(l - x) . (22) 

Here C = (x(t = 0)) =Lxp(x,t = 0)dx is the initial population state 
(McKane and Waxman, |2007 ). We plot $ and p(x, oo) in Figure 1 (red). 
In this case, we say that the peak-transition will never happen. The escape 
time is r = +oo. 

If, instead, there is additional factors such as (weak) mutation, the infinite 
escape time may not be a good estimation. Biologically, this mutation-like 
factor will constantly pull a population away from the monomorphic fixation 
state. It makes the substitution (peak-shift on the landscape) of a mutant 



possible, and the probability of such substitutions was often studied (Kimura 



1962). Mathematically, Eq. (20) causes an unexpected property of the escape 



time, that r/7i would change discontinuously with ANv (from +oo to 1 as 
ANu — > 1). We try to obtain better estimation for r by simulating the 
Wright-Fisher model and calculating the MFPT in the next section. 
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4. Escape time in infinite potential 

4-1- Two timescales: distribution's view 

To study the stochastic dynamics of a population, one may instead study 
the evolution of infinite many identical populations. In a diffusion model 
with two potential peaks, the distribution of the populations is expected to 



undergo two distinct stages of evolution in two different timescales (Barton 



and Rouhani, 1987). In the first timescale, the probability densities converge 
to the local potential peaks. Local equilibria are established on different po- 
tential peaks in 7i- In the second timescale, the probability densities transit 
(escape) between different equilibria. The escape rate of the probability den- 



sity in an attractive basin is assumed to obey an exponential law (Hanggi 



et al., 1990): 



Z (t) - Z (+oo) = T exp(-At) . (23) 

A is the average leaking rate of probability density. Here we define the cu- 
mulative probability density in the attractive basin (0, a) as 

Z (t)= [ p{x,t)dx. (24) 

Jo 

The subscript denotes the variables relating to the dynamical behaviors in 
(0,a). To is a time- independent function of system parameters. The leaking 
rate of Z (t) changes with time (actually, its change rate depends on the 
value of Z (t) at that time). We may characterize this process by the inverse 
of the average leaking rate of probability T2 = 1/A, which gives the timescale 
to establish the global equilibrium. 

Note that the leaking rate is the sum of contributions from both attractive 
basins (A = Ao + Ai = Tq 1 + r{~ ). Transitions in the two directions will 
eventually balance each other as t — > +00 at the global equilibrium: 

r -1 ■ Z (t = +00) = rf 1 ■Z l {t = +00) . (25) 

Here Zi(t) = 1 — Z (t) is the cumulative probability density in (a, 1). The 



escape time from the (0, a) basin is then given by (Hanggi et al. , 1990) 



T " = 1 - Ut = +oo) • (26) 

As the present landscape is a biological realization of the physical land- 
scape, we expect the biological population dynamics also obey this two-stage 
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evolution. To verify this, we simulate the change rate of P(t) (a discrete ver- 
sion of p(x, t)). The results in Figure [3] shows a clear two-timescale structure 
of landscape dynamics. With the probability densities initialized in (0, a), 
the decreasing rate of Z (t) is strictly exponential (Figure 111). Taking the log 
scale (in the sub-figure), the regressed slope of gives the escape rate toward 
the global equilibrium (A = A + Ai). The decay rate shows rigorous exponen- 
tial distribution, except for the sudden drop at the beginning period of time 
(~ 7i ~ 60 in the example, when the local equilibria has yet to be estab- 
lished). It validates the assumption of the two-scale dynamics, even though 
the potential peak at x = is singular and the established local equilibrium 
is not Gaussian. In the next section, to get an analytical estimation of r, we 
study the mean first passage time in infinite potential peak. 

4-2. MFPT in infinite potential 

In general diffusion model, one may calculate the mean first passage time 
from xq to Xi, by referring to the backward equation of Eq. (fl|. Without 
loss of generality, we study the stochastic jump out of the attractive basin 
(0,a). We study a population's MFPT through the valley point x = a to 
some state X\ > a, starting from xq ~ in (0, a). The interested interval is 
set as [0,xi], with x = the reflecting boundary and x = X\ the absorbing 



boundary (Gardiner, 1985): 



7mfpt(x ->■ xi) = / n / x ex P [ - ^{y)]dy / ex P [§{z)]dz . (27) 

Jx eL> {y) Jo 

Here $ is the potential landscape in Eq. ([7]). There is no requirement on the 
configuration of $ (finite peak, etc.) when applying Eq. (27). However, the 



MFPT is conceptually different from the escape time (Hanggi et al. , 1990). 
We will show how the MFPT can be used to analytically approximate the 
escape time in the present model. 

In a finite-barrier diffusion process, the escape time was shown to be ap- 
proximated by the MFPT and takes the form of the Arrhenius exponential 



factor (Gardiner, 1985). Previous approximation methods are mainly estab- 
lished on the following two assumptions: (1) A "sharp" valley around x = a 
on the landscape; (2) Gaussian-like probability distribution around x = 0. 
However, these two assumptions fail in the present model, as the landscapes 
have "fat" valleys and singular peaks under strong genetic drift (Figure 1). 



The equilibria established near a potential are not Gaussian (Barton and 
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Rouhani, 1987). In this section, we take another way to analytically approx- 
imate the escape time from the MFPT, in accordance with the landscape 
configuration in the present model. Under mutation and drift: 

Tmfpt(* -► *i) = AN r y-^{l - y)- 4N »dy f z m ^\l - zf N ^dz . 

Jx JO 

(28) 
The integral term z 4Nu x near z = accounts for the infinity of potential 
peak (and possibly infinite escape time) in the present model. We note that 
as ANv, ANfi — » the main contribution of the above integral comes from the 
inner integral in a small interval [0, y] (y < a), that is, the incomplete Beta 
function B(y;ANi>,ANfi). Under the same limit, it is numerically shown to 
be approximated by y/ANv. Thus the whole integral is approximately of a 



scale 1/v. More formally, we expand the incomplete Beta function in Eq. (28 ) 
under 0<1 — x% < 1 — y < 1 — z < 1 (expand (1 — z) iNfl ~ l near z = 0): 



B(y;4Nis,4Nfi)= [ z mv ~ l (l - z) m ^ l dz 
Jo 




V 

z 







\ n=0 fc=l 



,ANv °° „.n+ANv « u _ AM II 

= l + V - 17 ^ f 29) 

ANv ^n + iNv*-*- k v ; 

n=l fc=l 

The convergence of the expansion is obvious given < y < x\ < 1. Substitute 



5(y; 4iW, AN/j) and expand (1 — y) 4Nfl in the outer integral of Eq. (28 ), we 
obtain 

, v xi-xo , 4iV/i^x™ +1 -x ra+1 f r /A;-l + 4iV / i\ 

Zmfpt(zo -> a?i) = + > — M L + 

v v *-^ n+1 - LJ - \ k I 

71=1 fc = 2 v 7 

1 N ^(n + l)(n + ANu)l} 2 \ k J 

(30) 

The expansion converges under v > 0, fi < 1/AN. For the two limiting cases: 

(1) v — > 0: The expansion of Eq. (J29| becomes invalid. The leading term 

of the expansion changes from y^NvjANv to lny, which becomes sensitive 

to xo near then. To ensure the convergence of Tmfpt(^o - >■ #i) as Xq — > 0, 
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we need v 7^ 0; this is the condition for the escape problem (from x = 0) to 
be finite (as we have discussed in Section 3.3). On the other hand, we always 
have Tmfpt(0 — > x\) — >• 00 as v — > 0. 

(2) \x — >■ 1/AN: The expansion of (1 — y~)~ iN ^ would not converge for 
Xi —¥ 1, as the resulted series would then become a divergent harmonic 
series. This is also illustrated by the vanishing bi-stability of the system 



when AN n = 1 (Figure 1, yellow). To ensure the convergence of Eq. (30) as 
X\ — > 1, we need jj, < 1/4/V. 

4-3. Escape time 

From the results above, we are able to calculate the MFPT between any 
two points Xo and x\ that satisfy < xo < X\ < 1. The relation between 
the MFPT and the escape time, however, is not direct. The escape time, 
defined as the inverse of the average probability rate through the potential 



valley x = a (Kramers, 1940), cannot usually be (intuitively) estimated by 
T MFPT (0 — > a) on general potential energy landscape. This is because a 
population may still have some probability to return back to Xq immediately 
after its first arrival to x = a. In cases with axisymmetric landscape (with 
axis x = a), it is demonstrated that Tmfpt(0 — > a) should be compensated 



by a factor of 2 when approximating the escape time (Hanggi et al. 1990). 



Under 4iW, 4iV/i ^C 1, we have by Eq. ( 12 ) that a fa 1/2 and the approx- 



imately axisymmetric landscape. The escape time r is (taking xq = 0): 
t ^ 2 x T M fpt(0 -> a) 



r^w< 



v v ^— ' n + 1 J - J - V k 

n=\ k=2 



«d - «") e ,„ x1 r^ w „ n ^^ ■ (3D 



t;(»+i)(»+») 



fc=2 



Under ANv,AN[i ^C 1, the escape time is approximately independent of the 



initial state x in Eq. (31). The escape time is dominated by the leading 



term of the series 1/V, added by a remaining term of order AN \ijv: 

r w iz-^l + 1.23iV/i) , (32) 



The coefficient 1.23 is an approximation of the remaining series in Eq. (pi 
Under ANv,ANfi ^C 1, tq is thus much bigger than the relaxation time (< 
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2N) given in Eq. (19). This shows the separation of the two timescales and 
completes our inquiry for the escape time from (0, a). 

Another way to look at the MFPT in Eq. ( p0| is to set X\ — 1 and obtain 
the substitution time of A\ alleles. It differs from the escape time above by 
taking into account the dynamical details in the other attractive basin (a, 1): 

wo ., d = i + ^ ± 4- n ( k -^p^) + 

v v *-^ n + 1 - LJ - \ k I 

n=l fc=2 v 7 

^-^)E (n+1)( : +w 4 (^^)- 

(33) 

The necessary condition for its convergence (z/ > 0, fi < 1/AN) has been 
discussed in Section 3.1. In Appendix A we show that the condition is also 
sufficient. Biologically, we expect that Tmfpt(0 — > 1) > r , as a population 
would not have to arrive at x — 1 before it is identified to have escaped from 



(0, a). This difference is also mathematically demonstrated by Eqs. (31 )(33 ). 
Under the limit 4iW, 4Nfi <C 1, the two equations arrive at the same result 
7mfpt(0 — y 1) ~ r ~ 1/v. Numerical comparison of 2Tmfpt(0 — > a), 
7mfpt(0 — > a), our analytical approximation, and escape time simulated 
from the discrete model are given in Figure [5j 

In Figure [5j the escape time is best approximated by 2T M fpt(0 — > a). 
T MFPT (0 — y 1) also approximates the results well, but is always bigger. The 
estimation Tmfpt(0 — > a) is obviously not a good approximation: This is 
different from the usual model with finite sharp valley, where the escape 
time is approximated by T(0 — > X\) with x\ ~ a. The simulated escape time 
is always between T MFPT (0 — y a) and T M fpt(0 — > 1). Another observation 
is that 2T MF pt(0 — > a) is bigger/smaller than the simulated escape time 
when 4iW is small/big. The main reason is that the landscape is no longer 
symmetric when 4iW ^ 4Nfi (see Section 4.2 for more discussions). 

4-4- Models with weak selection 

With the effects of selection, the evolutionary rate is usually not linear- 
dependent on the gene frequency of A\. Its adaptive nature (if works alone) 
will drive a population monotonously to a fitness peak on the classical fitness 
landscape. Under mutation and genetic drift, a population is not expected 
to evolve towards a fitness peak, but towards a potential peak (Eq. ([l~6|)). 
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We study how the non-linear adaptive selection interacts with mutation and 
drift on the present potential landscape. 

Though f(x) no longer takes the linear form, the first timescale can still be 
estimated by 7i ~ 2N-0(1) under weak selection (ANs <C 1). This constraint 
makes sense in that most selections are weak as observed in nature. The 
general equation for the MFPT when there is mutation, drift and selection 



is obtained by substituting Eq.(lO) into Eq.(27) 



7mfpt(zo -> X\) 



AN 



XI 



:i-v)' 



ANfj, -ANv 



-27V 



[u(y)] dy 



■CO 



A-z 



,4JV>-1 ANv-l 



\U(Z 



2/V 



dz 



(34) 

The inner integral is no longer the standard incomplete Beta function. If we 
can expand the average fitness term [u(y)] 2N near the singular state x = 0, 
an analytical approximation for the MFPT can be obtained by combining 



the results with Eq. (33) 



An example of selection is the symmetric selection (Barton and Rouhani 



1987). The fitness setting is AiA t : A X A 2 : A 2 A 2 = 1 : 1— s : 1. Assuming the 



Hardy- Weinberg equilibrium, the effect of selective pressure can be described 



by the allele frequency (Gillespie, 1998). The average rate of change in x per 
generation by selection alone is M s (x) = —sx(l — x)(l — 2x). Here s is 
the fitness deficit of the heterozygote relative to the homozygote. We have 
1 — 2sx + 2sx 2 .li further assuming // = v, the potential landscape is 



</.' 



exactly axisymmetric to the axis a = 1/2 (also the potential valley): 
$(x) = (4Nfj, - 1) lnx(l - x) - ANsx + ANsx 2 , 



(35) 



plotted in Figure 1 (cyan). We study the effects of weak selection on the 
escape time on the basis of previous discussion. By expanding e - 4Nsx + 4Nsx 
near x = 0, the escape rate A 



r 1 is obtained as 



A = 1/(2 x T M fpt(0 -> a)) 



ANsy-4Nsy 2 (i 



y) 



AN l , 
'o 

/i/(l + 1.23Nfi + 0.67Ns) 



ry 

ANfi -4Nfii / e -4Nsz+4Nsz 2 M_ z \4Nn-l z 4N i i-l^ z 

Jo 

(36) 



Here a = 1/2 is the saddle (valley) state of the landscape in Eq. (35). We fix 



other parameters, give numerical comparisons among Eq. (34) (specified by 



Eq. (35) and take the inverse), Eq. (36) and discrete results in Figure [7j It 
can be noticed that our approximation also works for 4A^s ~ 1. 
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Another typical model considers the asymmetric selection and mutation. 
Under this condition, the fixation time of a beneficial or a deleterious muta- 
tion is often studied. We check the consistency of our results to the previous 
conclusion, and show that our results can be applied to more general cases 
with reverse mutation. For example, if we take s as the selective advantage 
oi A 1 over A 2 (s <C 1), such that M s = sx(l — x) (Kimura, 1964), the average 
fitness is given by w = 1 + 2sx . We study the average time of substitution 
of gene A\ by gene A 2 , which can be approximated by Tmfpt(0 — > 1). Note 
that to maintain a bi-stable system, we set 1/AN > /i, v. To take the expan- 



sion we further assume AN s < 1. Substitute above settings into Eq. (34) and 
obtain: 



T M fpt(0 ->> 1) « AN I 1+ UN(j,-4Ns\y 



1 1 - ANfi + ANs 



ANv 



+ 



l + 4iV^ 



dy 



(1 + 2Nfi - 2Ns) 
v 



(37) 



the substitution time of A x alleles. From this result, the selective advantage 
s decreases the substitution time approximately on a linear scale if ANs < 
1, consistent to the rate of substitution calculated under the same settings 
without backward mutations (/i = 0,s <C I, ANs < 1) (Gillespie, 1998): 



k 



1-e 



-2s 



~4Ni 



x 2Nv 



2Ns 



just the inverse of Eq.(37) if we take // = 0. If there is reverse mutation 



(// 7^ 0), the actual fixation would not happen. Our Eq. (37) can still be 
applied in such cases, carrying the meaning of general transition rate in the 
direction A\ — > A 2 . 



5. DISCUSSION 

5.1. Comparisons with previous work 

In the present work, we studied the two-timescale landscape dynamics 
in the Wright-Fisher diffusion model and estimated the escape time from 
the MFPT. In the 1-d diffusion model, actually, the time-dependent solu- 



tion p(x,t) of the diffusion equation can be analytically obtained (Kimura 



1964), and so is the explicit form of escape rate dZ (t)/dt. However, the 
exact analytical solution is in very complex form, and the structure of the 
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two-timescale dynamics is very difficult to identify from the complex math- 
ematics. Our landscape-based results, instead, enabled a very clear view of 
the separation of the two-stage evolution. The typical exponential decay of 



probability distribution in physics (Hanggi et al. , 1990) was also observed in 
this population genetics model. We employed the MFPT method to get an 
analytical approximation of the escape time To in the diffusion model, which 
is also found to be consistent to the eigenvalue estimations in the discrete 



Wright-Fisher model (Blythe and McKane, 2007). We overcame a technical 



difficulty for applying the classical method into the non-Gaussian equilibria 
with weak mutations, and further applied it in models with weak selections. 



From Eq. (32), the transition time is approximately independent of the 
population size 2N . It reminds us of Kimura's known rate formula for the 



neutral evolution (Kimura, 1962), or the rate of neutral substitution: 2Nv x 
1/2N = v, just the inverse of Eq. (32) if we take fi = 0. The inverse of the 



substitution rate gives the expected time of the appearance of a mutation 
that is destined to be fixed (Felsenstein, 2011). The coincidence of the two 



results happens under the limit 4iW <C 1, where the time for the desired 
mutant to be actually fixed (~ 2N) is negligible. For comparable v and 
1/2N, the population size iV will have significant effects on the transition 



rate. Our result Eq. (31) shows how the population size will have effect on 



the escape time. Moreover, Eq. (31) can further be applied under two-way 



mutations, which makes the fixation probability of a new mutant (and thus 
the rate of substitution) incalculable. 

Our results show that Kramers' classical escape time results derived from 
the MFPT or the flux-over-population method can be extended to the non- 
Gaussian distribution cases. Under 4iW, ANfi <C 1, the result does not show 
exponential dependency on the valley depth (or barrier height), but rather 
is controlled by the sharpness of the potential peak (see the sensitivities of 



Eq. (11) and Eq. (32) with respect to v). The difference is essentially induced 
by the special types of noise induced by genetic drift. On the other hand, 
under 4iW, AN [i C 1, we have 7i ~ 2N and J~2 ~ ^ _1 ; there is still the 



separation of different timescales, a natural result of our Eq. (31). 



The eigenvalue method was used to study the second example in our 
Section 4.4. The method failed to approximate the transition rate under very 
weak selection (s < 4/z), however, as the approach requires two "deterministic 
equilibria" (deterministic equilibrium is defined as a state which would be an 
equilibrium in the absence of random perturbations; here it can be simply 



considered as a mutation-selection landscape) (Barton and Rouhani, 1987). 
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However, the present Wright-Fisher model has a special form of noise of 
genetic drift (shown in Eq. (J3J)). When the noise is strong, it will result 
in two "stochastic equilibria" even when there is only one "deterministic 
equilibrium" . This is where the eigenvalue method fails. 

5.2. MFPT and escape time 

In general, an escaped population is assumed to be caught stable in an- 
other attractive basin (here (a, 1)). We do not take into account the cases 
in which a population comes back immediately after it escapes from the at- 
tractive basin (0, a). In models with a narrow potential barrier at x = a, 
the escape time was shown to be approximated by T MFPT (0 — > xi) with 



x\ > a (Gardiner, 1985). Here x\ is near a. However, this is not the case in 
the present model, as the potential valley is not narrow. Even after passing 
through some point a < X\ < 1, there may still be considerable possibil- 
ity for it to return to (0, a). In Figure [6j the dependence of the MFPT on 
the end point x = X\ is approximately linear. This conclusion requires that 



AN is, AN /i <C 1: From Eq. (30), the MFPT is dominated by the first linear 
term of x± — xq under this condition. From this near-linearity, we can see 
that the MFPT is conceptually very different from the real-time probability 



flow through each state (Hanggi et al. , 1990[ ), which should not be a constant 



value. For example, the probability flow (in the direction — y 1) near x = a 
should be much bigger than that near x = 0. Comparing with the escape 
time, the MFPT describes a transient event in the evolution, which has no 
direct implication for the future dynamics. An example is that, after first 
reaching x — a, there is still ~ 1/2 probability to go into (0, a) or (a, 1) then 
(see Appendix B for more detailed discussions). 

Figure [6] also numerically validates the relation r = 2 x T M fpt(0 — > a). 
The saddle x = a is a critical state in the population evolution: the most 
improbable state at long time observation and the unstable fixed point on 
the potential landscape. In a typical process with Gaussian local equilibria 
in physics and chemistry, the main escape time is spent in climbing over the 
saddle point. In the present Wright-Fisher model, the special form of the 
stochastic term is induced by genetic drift (Eq. pi), and the landscape config- 
uration is also different from those in the usual cases. The main difficulty of 
escape lies in overcoming the sharp peak at x — 0, where most probability 
densities are concentrated. 

Another information that can be read from Figure [6Jis that Tmfpt(0 — > 1) 
is generally bigger than the escape time. This is also easily observed from 
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the comparison between Eq. (31) and Eq. (33). This is because the escape 



time from (0, a) does not take into account the actual dynamics in the basin 
(a, 1) (Hanggi et al. , 1990), and a population does not have to reach x — 1 



before it is considered to have escaped. Eq. (31) is a better approximation 



5. 3. More on our potential landscape 

By theoretical analysis and computer simulations, we show that an evolv- 
ing population's probability distribution will first follow the uphill direction 
to the nearest potential peak in 7i, and then leak downhill through the sad- 
dle (valley) state and establish the global equilibrium in 7i- It verifies that 
the population dynamics are faithfully described by the present landscape. 
Its coherency is also demonstrated in the limiting case of pure genetic drift 
as shown in Section 3.3 and Figure 1 (Green) with Z = +oo. As shown 
in Section 2.1, $ relates to p(x, t = +oo) through the Boltzmann-Gibbs 
distribution (if Z < +oo), but is essentially a dynamical description of the 
system. Its validity does not require a normalizable equilibrium distribution. 

The potential landscape Eq. d7|) can be compared to the classical fitness 
landscape, which presents only the effects of selection. Other biological fac- 
tors may generate various evolutionary mechanisms on the fitness landscape 



without a unified description, along with other controversies (Kaplan, 2008). 
Also, by only taking the measure of fitness, there may be inconsistencies 
between the dynamics and biology. Under the effects of other forces (e.g. 
mutation and drift), a population is not always expected to move toward the 
fitness peak, but toward the potential peak (Section 3.2). One example is the 
term "neutral evolution" commonly used in the absence of selection, where 
different allele-frequency states of a population are not necessarily equally 
favored by evolution (except the special case 4Nu = ANfi = 1), shown 
in Figure 1. The present potential landscape (as also noticed in literature 



(Burger 2000)) may serve as a substitute for Wright's original landscape that 
visualizes and quantifies the evolutionary process in a globally coherent way. 
An extension to the fitness landscape is the "deterministic landscape" 
(Barton and Rouhani, 1987), which integrates all deterministic factors of 
evolution but does not consider genetic drift. The uphill rate on the de- 
terministic landscape is given by M(x) (instead of f(x), compared to our 
potential landscape). For example, under mutation and random drift, we 
have M(x) = —fix + z/(l — x) . It has a zero point om = v /{ v + aO? which 
is also the unique peak state of the deterministic landscape. It fails to de- 
scribe the Wright-Fisher model when genetic drift is strong (4 Nv, 4Np < 1 , 



our potential landscape is bi-peaked); the associated approaches also fail for 
such cases (see Sections 4.1 and 3.3). Our directed rate f(x) in Eq. (13), 
instead, gives the correct expectation of directed evolution. The key point 
here is that genetic drift (noise) also contributes to the directed force in the 
long-term observations. 

Another extension of the classical fitness landscape is the free fitness func- 



tion (Barton and Coe, 2009) in consideration of the analogy with thermo- 



dynamics. It uses the normalization constant Z as a generating function of 
system's macroscopic information and requires normalizable probability dis- 
tribution. Moreover, the associated maximum entropy approximation fails 
under weak mutations. Our present framework does not have certain con- 
straints, and the validity of our landscape construction and the associated 
approaches is tested in the whole relevant parameter space. It has already 



been applied in the study of Muller's ratchet (Jiao and Ao, 2012), a special 
case where no backward mutations exist. 



5-4- Normalization constants and fixation 
By taking v 



in Eq. (30), we have T\ 



foo. No escape is expected 
to happen once a population "trapped" into the neighborhood of x = 0. 



In Eq. (27), the impossibility comes essentially from the infinity of the in- 
complete Beta function B(y; 4Ni>, ANfi) in Eq. (29). More formally, if we 



define a partial normalization constant for each attractive basin (taking the 
mutation-drift case as an example) as 



x 



ANv-l, 



1 _ x ) iN ^dx 



Zx 



ANv-X, 



1 



) 4N »- l dx 



(38) 
(39) 



Note that Zq can also be obtained from the definition in Eq. (24) by taking 
t = +oo: Zq = Zo(t = +oo). Here Z\ = 1 — Zq is similar. The mathematical 
condition for the biological fixation at x = (or x = 1) should be Z = +oo 
(or Zi = +oo). When combined with previous discussions, we conclude that 
$(0) = +oo (or $(1) = +oo) does not necessarily imply fixation of A\ (or 
Aq) genes. The condition for a population starting from any initial state 
to be fixed at a monomorphic gene state x = or x = 1 is determined by 
Z = +oo, Z 1 < +oo or Z < +oo, Z\ = +oo. If Z = Z x = +oo, the 
fixation will happen at either x = or x = 1 on chance. Another obser- 
vation from the present results is the emerging of absorbing boundaries at 
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the fixation state; the boundary conditions "artifically" set by (McKane and 



Waxman, 2007) are more naturally and generally derived then. Note that 



certain boundary conditions have also been discussed using MFPT in the 
literature of stochastic theory (Karlin and Taylor, 1981). However, without 



a proper framework or detailed discussions on the use of MFPT, the results 
may not be rigorous and may even be erroneous in some cases. A key issue 
here is the relation between the escape time and MFPT: the MFPT from to 
1/2 does not necessary mean the escape from the (0, a) (|Karlin and Taylor 



1981), as shown in our Figure [5] and [7J Our last comment is that unnor- 
malizable distributions in the diffusion model do not generate real problems 
for understanding the original discrete model. It instead provides important 
dynamical and equilibrium information for the understanding of the system. 
We summarize above conclusions in Table 1. 



Table 1: Summary of the observations in Section 4.3. Zq and Z\ are the partial normal- 
ization constants defined in Eqs. (38) and (39). r and t\ are the respective escape times. 



The "Absorb-bound." column gives where the absorbing boundary emerges. 
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< oo 


x = 




x = 
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x = 1 




x = 1 
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= oo 


= oo 


= oo 


x = or x = 


1 


a; = 0,1 



5.5. Comments on the "stochastic tunneling" 

In study of a three-phase transition problem, a "stochastic tunneling" 
effect was termed that allows transition from one state to another, without 



passing through the middle state (Iwasa et al. , 2004). In light of the present 



framework of potential landscape and discussions of escape time, we com- 
mented that "stochastic tunneling" is a misused term. It is clear that the 
essential feature of the tunneling effect does not take place here: Tunneling 
is a quantum mechanical effect. It is tied to the laws of wave mechanics 



going through under the potential barrier (Anslyn and Dougherty, 2006 Ao 



and Thouless, 1994). A potential barrier (potential valley) is not to be over- 
come but is tunneled through, which is classically impossible. Furthermore, 
the tunneling effect is approximately temperature- independent, while the 
"stochastic tunneling" disappears when temperature (noise intensity, or here 
the population size) decreases to zero (population size increases to infinity). 
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The first step of fixation of the deleterious mutation would never happen 
with zero noise. In fact, this process is just an "old-type" saddle-passage es- 
cape event on an potential landscape. The proper term should be "thermal 
activation" . 
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Appendix A. Convergence of Eq. (33) 



sum 



Under v > 0, the convergence of Eq. (33) relies on the convergence of the 

fc-l + 4N>\ 1 



s=£n 



k 



n + l 



n=2 k=2 

We use Raabe's test for series convergence from standard textbooks of real 
analysis. For < ANfi < 1, we denote 



n 

fc=2 



jfe-l 



A- 



ANu\ 



J n + l 



Obviously c n is positive for all n > 0. First, we have 

c n+l 



lim 

n-K» c r 

We then calculate the Raabe terms 



1 . 



(A.l) 



i^ = n(^-l) = (4iV>-2) 



n 



n + 2 



Here 4iV/i — 2 is a constant less than — 1. By taking the limit n — > oo, 

lim R n = 4Nfi - 2 < -1 (A.2) 



The two conclusions in Eqs.(A.l, A.2) verify the convergence of the partial 
sum S n under < 4Nft < 1. 

Appendix B. Interpretation of the factor 2 

The choice of factor 2 is because a population will have (on average) 1/2 
probability to be caught stable in the (a, 1) basin then. This is not always the 
truth, though, as valley X\ = a is chosen as a perfect absorbing boundary 



(sink) rather than a smooth distribution of sinks in (a, 1) (Hanggi et al 



1990). For the factor of 2 to be exact, we need the valley point to have 



zero slope ($(s = a)) and the landscape to be axisymmetric near the valley. 
Asymmetry of the landscape far from the valley state will bring higher-order 
errors to the factor of 2, which is neglected for the present concern and needs 
further investigations. On the other hand, differences between 2T mfpt and 
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the simulated escape time may also come from the diffusion approximation 
and expansion of the MFPT. 

One detailed interpretation of this factor 2 in a limiting case is given 
below. Assume that the rates of uphill (~ 7i) and downhill (T MF pt(0 —> a), 
denoted as T a ) movements are well separated (e.g. 4Nv,4:Nfi <C 1). Once 
a population reaches x = a in T a , it has (on average) probability 1/2 to fall 
into either attractive basin and immediately (~ 7i, assumed to be negligible 
compared to T a ) reaches the corresponding potential peak; once falling back 
to x = 0, it will wait another T a to reach x = a again; it then again has 1/2 
chance to reach x = 1 or return to x = immediately. Assume this process 
continues. The expected time to leave x = can then be obtained as 

r^T a xi + 2T a xl + ...+nT a xl + ..., 
= 2T a , (B.l) 



the same result as Eq. (31) 



Appendix C. Discrete Wright-Fisher process 

The original Wright-Fisher model is discrete both in time (number of 
generations) and space (number of copies of Ai). It considers the evolution 
of the probability distribution function P t (a vector of 2 A + 1 elements) with 
time t. The t th generation sampled 2N times to give the t+ 1 th generations. 
The probability that these 2A trials of sampling a population with i copies 
of Ai gene will give j copies of A\ gene is given by the (i, j)th element of the 
transition probability matrix M, defined by the binomial distribution: 

M tJ =C] N p(^(l-p( l )) 2N ^ . (C.l) 

C 2N is the number of combinations to choose j genes from a gene pool of 
size 2A^. p(i) is the probability of choosing a A\ gene from the pool. To 
give the explicit form of p(i), we denote y = i/2N. p(i) is determined by the 
biological factors considered in the model. Here under genetic drift, mutation 
and selection: 

(WW + W 12 y(l - y))(l -fi) + (Ml -y) + W 22 (l - yf)u 

P[l) ~ WW + 2W 12 y(l -y) + W 22 {1 - y) 2 

(C.2) 
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is the probability of sampling an Ai gene in the population gene pool. The 
system evolves as 

P t+1 = P t -M (C.3) 

In the simulation in Section 4.1, the model is under mutation and genetic 
drift, and the fitnesses are chosen as Wu : W\ 2 : W 22 = 1:1:1. In the first 
example of Section 4.4, the fitnesses are Wu : Wi 2 : W 2 2 = 1 : 1 — s : 1. 
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Figure 1: Potential landscapes and corresponding equilibrium distributions under different 
parameter settings in the Wright-Fisher model, differentiated by both the colors and Ro- 
man indexes. In all cases there is N — 50. The following five colored landscape contours 
are generated from Eq. (11) under mutation and genetic drift: Red (I): /j, = v = 0. 
(II): fi = 0.0005, v = 0.001. Blue (III): /z = v = 0.005. Yellow (IV): fi = 0.005, v = 
Magenta (V): fj, = 0.01, v — 0.001. The last one is generated from Eq. (35), considering 
mutation, drift, and selection: Cyan (VI): fi = v = 0.002, s = 0.1. The two red arrows in 
(b) denote the Dirac delta functions. 27 
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Figure 2: Visualization of the two-scale dynamics on a typical bi-stable landscape in the 
present model. It gives the most probable state of a population (denoted as a balloon, 
which always searches for a higher "altitude" to stay) in different timcscalcs visualized on 
the potential landscape. The parameters are: 2N = 20, jjl = 0.0005, v = 0.0015. 
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Figure 3: Simulations of the discrete Wright-Fisher model under mutation and random 
drift, x-axis gives the number of A\ alleles and y-axis is the probability distribution. 
Parameter settings: 2N = 20, fx = 0.0005, v = 0.0015, so that T x « 20, r sw 666. (a) 
shows that the initial state is set to x — 0.2. (e) and (f) show the establishment of the 
equilibrium distribution after long enough time. 
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Figure 4: Simulation of the escape rate from the attractive basin (0, a) under mutation 
and drift. Parameter settings: N = 30, // = 0.0001,^ = 0.0005. The main plot de- 
scribes how the cumulative probability density in (0,a) attractive basin Zo(t) (defined 
in Eq. (24)) changes with time. The inserted figure is the value —\n[Z (t) — Z (+oo)], 
whose slope gives the flux rate between the two attractive basins. The zoomed-in fig- 
ure shows the same value in the first timescale. The simulated values are: T\ sa 
62. 77(regressing the steady exponential interval (0, JV)), T 2 sw 1464, f « 1755. Under 
the same setting, the theoretical expectations are: T\ = 60, T 2 — 1666, t = 2000) 
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Figure 5: Comparison between two times the MFPT from to a (solid), our analytical 
approximation (dashed), the MFPT from to 1 (point), the MFPT from to a (dot- 
ted), and the simulated escape time (crosses) of the discrete Wright-Fisher model under 
mutation and genetic drift. Parameter settings: N = 100, /i = 0.00005. 
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Figure 6: Comparisons between the MFPT from to xi (solid), our analytical approxima- 
tion of the MFPT (dashed), and the discrete MFPT calculated from the Master equation 
(square). The parameter setting is N = 100, v = 0.00025, \x = 0.00001. The saddle point 
a = 0.4747. The crosses are the simulated escape time of the discrete Wright-Fisher model 
under the same setting. 
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Figure 7: Comparison between two times the MFPT from to a (solid), our analytical ap- 
proximation (dashed), the MFPT from to 1 (point), the MFPT from to a (dotted), and 
the simulated escape time (crosses) of the discrete Wright-Fisher model under selection, 
mutation, and genetic drift. Parameter settings: N — 50,^ = 0.00025, \i = 0.00001. 
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